!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!                                                             
!   glint_main.F90 - part of the Community Ice Sheet Model (CISM)  
!                                                              
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!   Copyright (C) 2005-2018
!   CISM contributors - see AUTHORS file for list of contributors
!
!   This file is part of CISM.
!
!   CISM is free software: you can redistribute it and/or modify it
!   under the terms of the Lesser GNU General Public License as published
!   by the Free Software Foundation, either version 3 of the License, or
!   (at your option) any later version.
!
!   CISM is distributed in the hope that it will be useful,
!   but WITHOUT ANY WARRANTY; without even the implied warranty of
!   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!   Lesser GNU General Public License for more details.
!
!   You should have received a copy of the Lesser GNU General Public License
!   along with CISM. If not, see <http://www.gnu.org/licenses/>.
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#ifdef HAVE_CONFIG_H
#include "config.inc"
#endif

module glint_main

  !>  This is the main glimmer module, which contains the top-level 
  !>  subroutines and derived types comprising the glimmer ice model.

  use glimmer_global, only: dp, fname_length
  use glint_type
  use glint_global_grid
  use glad_constants
  use glint_anomcouple

  use glimmer_paramets, only: stdout, GLC_DEBUG

  implicit none

  ! ------------------------------------------------------------
  ! glint_params derived type definition
  ! This is where default values are set.
  ! ------------------------------------------------------------

  type glint_params 

     !> Derived type containing parameters relevant to all instances of 
     !> the model - i.e. those parameters which pertain to the global model. 

     ! Global grids used ----------------------------------------

     type(global_grid) :: g_grid      !> The main global grid, used for 
                                      !> input and most outputs
     type(global_grid) :: g_grid_orog !> Global grid used for orography output.

     ! Ice model instances --------------------------------------

     integer                                   :: ninstances = 1       !> Number of ice model instances
     type(glint_instance),pointer,dimension(:) :: instances  => null() !> Array of glimmer\_instances

     ! Global model parameters ----------------------------------

     integer  :: tstep_mbal = 1        !> Mass-balance timestep (hours)
     integer  :: start_time            !> Time of first call to glint (hours)
     integer  :: time_step             !> Calling timestep of global model (hours)

     ! Parameters that can be set by the GCM calling Glint

     logical  :: gcm_smb = .false.     !> If true, receive surface mass balance from the GCM 
     logical  :: gcm_restart = .false. !> If true, restart the model from a GCM restart file
     character(fname_length) :: gcm_restart_file   !> Name of restart file
     integer  :: gcm_fileunit = 99     !> Fileunit specified by GCM for reading config files
   
     ! Averaging parameters -------------------------------------

     integer  :: av_start_time = 0   !> Holds the value of time from 
                                     !> the last occasion averaging was restarted (hours)
     integer  :: av_steps      = 0   !> Holds the number of times glimmer has 
                                     !> been called in current round of averaging.
     integer  :: next_av_start = 0   !> Time when we expect next averaging to start
     logical  :: new_av     = .true. !> Set to true if the next correct call starts a new averaging round

     ! Averaging arrays -----------------------------------------

     real(dp),pointer,dimension(:,:) :: g_av_precip  => null()  !> globally averaged precip
     real(dp),pointer,dimension(:,:) :: g_av_temp    => null()  !> globally averaged temperature 
     real(dp),pointer,dimension(:,:) :: g_max_temp   => null()  !> global maximum temperature
     real(dp),pointer,dimension(:,:) :: g_min_temp   => null()  !> global minimum temperature
     real(dp),pointer,dimension(:,:) :: g_temp_range => null()  !> global temperature range
     real(dp),pointer,dimension(:,:) :: g_av_zonwind => null()  !> globally averaged zonal wind 
     real(dp),pointer,dimension(:,:) :: g_av_merwind => null()  !> globally averaged meridional wind 
     real(dp),pointer,dimension(:,:) :: g_av_humid   => null()  !> globally averaged humidity (%)
     real(dp),pointer,dimension(:,:) :: g_av_lwdown  => null()  !> globally averaged downwelling longwave (W/m$^2$)
     real(dp),pointer,dimension(:,:) :: g_av_swdown  => null()  !> globally averaged downwelling shortwave (W/m$^2$)
     real(dp),pointer,dimension(:,:) :: g_av_airpress => null() !> globally averaged surface air pressure (Pa)
     real(dp),pointer,dimension(:,:,:) :: g_av_qsmb => null()   ! globally averaged surface mass balance (kg m-2 s-1)
     real(dp),pointer,dimension(:,:,:) :: g_av_tsfc => null()   ! globally averaged surface temperature (deg C)
     real(dp),pointer,dimension(:,:,:) :: g_av_topo => null()   ! globally averaged surface elevation   (m)

     ! Fractional coverage information --------------------------
     ! Note: these are only valid on the main task
     real(dp),pointer,dimension(:,:) :: total_coverage  => null()     !> Fractional coverage by 
                                                                      !> all ice model instances.
     real(dp),pointer,dimension(:,:) :: total_cov_orog  => null()     !> Fractional coverage by 
                                                                      !> all ice model instances (orog).
     logical                         :: coverage_calculated = .false. !> Have we calculated the
                                                                      !> coverage map yet?
     ! File information -----------------------------------------

     character(fname_length) :: paramfile      !> Name of global parameter file

     ! Accumulation/averaging flags -----------------------------

     logical :: need_winds=.false. !> Set if we need the winds to be accumulated/downscaled
     logical :: enmabal=.false.    !> Set if we're using the energy balance mass balance model anywhere

     ! Anomaly coupling for global climate ------------------------------------------

     type(anomaly_coupling) :: anomaly_params !> Parameters for anomaly coupling

  end type glint_params

  ! Private names -----------------------------------------------

  private glint_allocate_arrays
  private glint_readconfig, calc_bounds, check_init_args
  private compute_ice_sheet_grid_mask
  private compute_icemask_coupled_fluxes

    !---------------------------------------------------------------------------------------
    ! Some notes on coupling to the Community Earth System Model (CESM).  These may be applicable
    ! for coupling to other GCMs:
    !
    ! When coupled to CESM, Glint receives three fields from the coupler on a global grid
    ! in each of several elevation classes:
    !   qsmb = surface mass balance (kg/m^2/s)
    !   tsfc = surface ground temperature (deg C)
    !   topo = surface elevation (m)
    ! Both qsmb and tsfc are computed in the CESM land model.
    ! Five fields are returned to CESM on the global grid.
    !   gfrac = fractional ice coverage
    !   gtopo = surface elevation (m)
    !   ghflx = heat flux from the ice interior to the surface (W/m^2)
    !   grofi = ice runoff (i.e., calving) (kg/m^2/s)
    !   grofl = liquid runoff (i.e., basal melting; the land model handles sfc runoff) (kg/m^2/s)
    ! The first three (gfrac, gtopo, ghflx) are returned for each elevation class of each grid cell;
    ! the last two are returned only for the full gridcell.
    ! The land model has the option to update its ice coverage and surface elevation, given
    ! the fields returned from Glint.
    !
    ! There are two driver subroutines in this module for CESM coupling: 
    !  initialise_glint_gcm (for initialization) and glint_gcm (for timestepping).
    ! These drivers loop over the ice sheet model instances (just Greenland for now,
    !  but will simulate Antarctica and other ice sheets later).
    !  
    ! The other driver subroutines, based on the original Glint code in Glimmer version 1,
    !  are initialise_glint and glint.
    ! These subroutines are usually run with temp (= 2-m air temperature) and precip as input.
    !  The surface mass balance is computed using a daily or annual PDD scheme.
    ! It would be possible to call these subroutines from CESM and use the PDD scheme,
    !  but this option has not been tested.
    !
    !---------------------------------------------------------------------------------------

contains

  !TODO - Try calling subroutine initialise_glint from CESM to estimate SMB in PDD mode?
  !       We would have only one elevation class per grid cell and would not return upscaled fields.

  subroutine initialise_glint(params,                         &
                              lats,         longs,            &
                              time_step,    paramfile,        &
                              latb,         lonb,             &
                              orog,         albedo,           &
                              ice_frac,     veg_frac,         &
                              snowice_frac, snowveg_frac,     &
                              snow_depth,                     &
                              orog_lats,    orog_longs,       &
                              orog_latb,    orog_lonb,        &
                              output_flag,  daysinyear,       &
                              snow_model,   ice_dt,           &
                              extraconfigs, start_time,       &
                              gmask,                          &
                              gcm_restart,  gcm_restart_file, &
                              gcm_debug,    gcm_fileunit)

    !> Initialises the model
    !> For a multi-processor run, the main task should specify lats & longs spanning
    !> the full global domain; the other tasks should give 0-size lats & longs arrays
    !> Output arrays on the global grid are only valid on the main task

    use glimmer_config
    use glint_initialise
    use glimmer_log
    use glimmer_filenames
    use glint_upscale, only: glint_upscaling
    use parallel, only: main_task
    implicit none

    ! Subroutine argument declarations --------------------------------------------------------

    type(glint_params),              intent(inout) :: params      !> parameters to be set
    real(dp),dimension(:),           intent(in)    :: lats,longs  !> location of gridpoints 
                                                                  !> in global data.
    integer,                         intent(in)    :: time_step   !> Timestep of calling model (hours)
    character(*),dimension(:),       intent(in)    :: paramfile   !> array of configuration filenames.
    real(dp),dimension(:),  optional,intent(in)    :: latb        !> Locations of the latitudinal 
                                                                  !> boundaries of the grid-boxes.
    real(dp),dimension(:),  optional,intent(in)    :: lonb        !> Locations of the longitudinal
                                                                  !> boundaries of the grid-boxes.
    real(dp),dimension(:,:),optional,intent(out)   :: orog        !> Initial global orography
    real(dp),dimension(:,:),optional,intent(out)   :: albedo      !> Initial albedo
    real(dp),dimension(:,:),optional,intent(out)   :: ice_frac    !> Initial ice fraction 
    real(dp),dimension(:,:),optional,intent(out)   :: veg_frac    !> Initial veg fraction
    real(dp),dimension(:,:),optional,intent(out)   :: snowice_frac !> Initial snow-covered ice fraction
    real(dp),dimension(:,:),optional,intent(out)   :: snowveg_frac !> Initial snow-covered veg fraction
    real(dp),dimension(:,:),optional,intent(out)   :: snow_depth  !> Initial snow depth 
    real(dp),dimension(:),  optional,intent(in)    :: orog_lats   !> Latitudinal location of gridpoints 
                                                                  !> for global orography output.
    real(dp),dimension(:),  optional,intent(in)    :: orog_longs  !> Longitudinal location of gridpoints 
                                                                  !> for global orography output.
    real(dp),dimension(:),  optional,intent(in)    :: orog_latb   !> Locations of the latitudinal 
                                                                  !> boundaries of the grid-boxes (orography).
    real(dp),dimension(:),  optional,intent(in)    :: orog_lonb   !> Locations of the longitudinal
                                                                  !> boundaries of the grid-boxes (orography).
    logical,                optional,intent(out)   :: output_flag !> Flag to show output set (provided for
                                                                  !> consistency)
    integer,                optional,intent(in)    :: daysinyear  !> Number of days in the year
    logical,                optional,intent(out)   :: snow_model  !> Set if the mass-balance scheme has a snow-depth model
    integer,                optional,intent(out)   :: ice_dt      !> Ice dynamics time-step in hours
    type(ConfigData),dimension(:),optional ::  extraconfigs !> Additional configuration information - overwrites
                                                                  !> config data read from files
    integer,                optional,intent(in)    :: start_time  !> Time of first call to glint (hours)
    integer, dimension(:,:),  optional,intent(in)  :: gmask       !> mask = 1 where global data are valid
    logical,                  optional,intent(in)  :: gcm_restart ! logical flag to restart from a GCM restart file
    character(*),             optional, intent(in) :: gcm_restart_file ! restart filename for a GCM restart
                                                                  ! (currently assumed to be CESM)
    logical,                  optional,intent(in)  :: gcm_debug   ! logical flag from GCM to output debug information
    integer,                  optional,intent(in)  :: gcm_fileunit! fileunit for reading config files

    ! Internal variables -----------------------------------------------------------------------

    type(ConfigSection), pointer :: global_config, instance_config, section  ! configuration stuff
    character(len=100) :: message                                            ! For log-writing
    character(fname_length),dimension(:),pointer :: config_fnames=>null()    ! array of config filenames

    integer :: i, j
    real(dp),dimension(:,:),allocatable :: orog_temp, if_temp, vf_temp, sif_temp,  &
                                           svf_temp,  sd_temp, alb_temp      ! Temporary output arrays
    integer,dimension(:),allocatable :: mbts,idts ! Array of mass-balance and ice dynamics timesteps
    logical :: anomaly_check ! Set if we've already initialised anomaly coupling

    if (present(gcm_debug)) then
       GLC_DEBUG = gcm_debug
    endif

    if (GLC_DEBUG .and. main_task) then
       write(stdout,*) 'Starting initialise_glint'
    end if

    ! Initialise start time and calling model time-step ----------------------------------------
    ! We ignore t=0 by default 

    params%time_step = time_step

    if (present(start_time)) then
       params%start_time = start_time
    else
       params%start_time = time_step
    end if

    params%next_av_start = params%start_time

    ! Initialisation for runs coupled to a GCM
 
    params%gcm_restart = .false.
    if (present(gcm_restart)) then
       params%gcm_restart = gcm_restart
    endif

    params%gcm_restart_file = ''
    if (present(gcm_restart_file)) then
       params%gcm_restart_file = gcm_restart_file
    endif

    params%gcm_fileunit = 99
    if (present(gcm_fileunit)) then
       params%gcm_fileunit = gcm_fileunit
    endif

!    nec = 1
!    if (present(gcm_nec)) then
!       nec = gcm_nec
!    endif

    if (GLC_DEBUG .and. main_task) then
       write(stdout,*) 'time_step =', params%time_step
       write(stdout,*) 'start_time =', params%start_time
       write(stdout,*) 'next_av_start =', params%next_av_start
    end if

    ! Initialise year-length -------------------------------------------------------------------

    if (present(daysinyear)) then
       call glad_set_year_length(daysinyear)
    end if

    if (GLC_DEBUG .and. main_task) then
       write(stdout,*) 'Initialize global grid'
       write(stdout,*) 'present =', present(gmask)
    end if

    ! Initialise main global grid --------------------------------------------------------------

    if (present(gmask)) then
       call new_global_grid(params%g_grid, longs, lats, lonb=lonb, latb=latb, mask=gmask)
    else
       call new_global_grid(params%g_grid, longs, lats, lonb=lonb, latb=latb)
    endif

    if (GLC_DEBUG .and. main_task) then
       write (stdout,*) ' ' 
       write (stdout,*) 'time_step (hr)  =', params%time_step
       write (stdout,*) 'start_time (hr) =', params%start_time
       write (stdout,*) 'Called new_global_grid '
       write (stdout,*) 'g_grid%nx =', params%g_grid%nx
       write (stdout,*) 'g_grid%ny =', params%g_grid%ny
       write (stdout,*) ' '
       write (stdout,*) 'g_grid%lons =', params%g_grid%lons
       write (stdout,*) ' '
       write (stdout,*) 'g_grid%lats =', params%g_grid%lats
       write (stdout,*) ' '
       write (stdout,*) 'g_grid%lon_bound =', params%g_grid%lon_bound
       write (stdout,*) ' '
       write (stdout,*) 'g_grid%lat_bound =', params%g_grid%lat_bound
       do j = 5, 10
          write (stdout,*)
          write (stdout,*) 'j, g_grid%mask =', j, params%g_grid%mask(:,j)
       enddo
    end if

    ! Initialise orography grid ------------------------------------

    call check_init_args(orog_lats, orog_longs, orog_latb, orog_lonb)

    if (present(orog_lats) .and. present(orog_longs)) then
       call new_global_grid(params%g_grid_orog, orog_longs, orog_lats,  &
                            lonb=orog_lonb, latb=orog_latb)
    else
       call copy_global_grid(params%g_grid, params%g_grid_orog)
    end if

    ! Allocate arrays -----------------------------------------------

    call glint_allocate_arrays(params)

    ! Initialise arrays ---------------------------------------------

    params%g_av_precip  = 0.d0
    params%g_av_temp    = 0.d0
    params%g_max_temp   = -1000.d0
    params%g_min_temp   =  1000.d0
    params%g_temp_range = 0.d0
    params%g_av_zonwind = 0.d0
    params%g_av_merwind = 0.d0
    params%g_av_humid   = 0.d0
    params%g_av_lwdown  = 0.d0
    params%g_av_swdown  = 0.d0
    params%g_av_airpress= 0.d0

    ! ---------------------------------------------------------------
    ! Zero coverage maps and normalisation fields for main grid and
    ! orography grid
    ! ---------------------------------------------------------------

    params%total_coverage = 0.d0
    params%total_cov_orog = 0.d0

    if (GLC_DEBUG .and. main_task) then
       write(stdout,*) 'Read paramfile'
       write(stdout,*) 'paramfile =', paramfile
    end if

    ! ---------------------------------------------------------------
    ! Determine how many instances there are, according to what
    ! configuration files we've been provided with
    ! ---------------------------------------------------------------

    if (size(paramfile) == 1) then
       ! Load the configuration file into the linked list
       call ConfigRead(process_path(paramfile(1)), global_config, params%gcm_fileunit)    
       call glint_readconfig(global_config, params%ninstances, config_fnames, paramfile) ! Parse the list
    else
       params%ninstances = size(paramfile)
       allocate(config_fnames(params%ninstances))
       config_fnames = paramfile
    end if

    allocate(params%instances(params%ninstances))
    allocate(mbts(params%ninstances), idts(params%ninstances))

    if (GLC_DEBUG .and. main_task) then
       write(stdout,*) 'Number of instances =', params%ninstances
       write(stdout,*) 'Read config files and initialize each instance'
    end if

    ! ---------------------------------------------------------------
    ! Read config files, and initialise instances accordingly
    ! ---------------------------------------------------------------

    call write_log('Reading instance configurations')
    call write_log('-------------------------------')

    anomaly_check = .false.

    do i = 1, params%ninstances

       call ConfigRead(process_path(config_fnames(i)),instance_config, params%gcm_fileunit)
       if (present(extraconfigs)) then
          if (size(extraconfigs)>=i) then
             call ConfigCombine(instance_config,extraconfigs(i))
          end if
       end if

       call glint_i_initialise(instance_config,    params%instances(i),     &
                               params%g_grid,      params%g_grid_orog,      &
                               mbts(i),            idts(i),                 &
                               params%need_winds,  params%enmabal,          &
                               params%start_time,  params%time_step,        &
                               params%gcm_restart, params%gcm_restart_file, &
                               params%gcm_fileunit)

       params%total_coverage = params%total_coverage + params%instances(i)%frac_coverage
       params%total_cov_orog = params%total_cov_orog + params%instances(i)%frac_cov_orog

       ! Initialise anomaly coupling
       if (.not.anomaly_check) then 
          call anomaly_init(params%anomaly_params, instance_config)
          if (params%anomaly_params%enabled .and. &
               (params%anomaly_params%nx/=params%g_grid%nx .or. &
                params%anomaly_params%ny/=params%g_grid%ny) ) then
             call write_log("Anomaly coupling grids have different "// &
                  "sizes to GLINT coupling grids",GM_FATAL,__FILE__,__LINE__)
          end if
          if (params%anomaly_params%enabled) anomaly_check=.true.
       end if

    end do

    ! Check that all mass-balance time-steps are the same length and 
    ! assign that value to the top-level variable

    params%tstep_mbal = check_mbts(mbts)
    if (present(ice_dt)) then
       ice_dt = check_mbts(idts)
    end if

    if (GLC_DEBUG .and. main_task) then
       write(stdout,*) 'tstep_mbal =', params%tstep_mbal
       write(stdout,*) 'start_time =', params%start_time
       write(stdout,*) 'time_step =',  params%time_step
       if (present(ice_dt)) write(stdout,*) 'ice_dt =', ice_dt
    end if

    ! Check time-steps divide into one another appropriately.

    if (.not.(mod(params%tstep_mbal,params%time_step)==0)) then
       print*,params%tstep_mbal,params%time_step
       call write_log('The mass-balance timestep must be an integer multiple of the forcing time-step', &
            GM_FATAL,__FILE__,__LINE__)
    end if

    ! Check we don't have coverage greater than one at any point.

    where (params%total_coverage > 1.d0) params%total_coverage = 1.d0
    where (params%total_cov_orog > 1.d0) params%total_cov_orog = 1.d0
    params%coverage_calculated=.true.

    ! Zero optional outputs, if present

    if (present(orog))         orog = 0.d0
    if (present(albedo))       albedo = 0.d0
    if (present(ice_frac))     ice_frac = 0.d0
    if (present(veg_frac))     veg_frac = 0.d0
    if (present(snowice_frac)) snowice_frac = 0.d0
    if (present(snowveg_frac)) snowveg_frac = 0.d0
    if (present(snow_depth))   snow_depth = 0.d0

    ! Allocate arrays

    allocate(orog_temp(params%g_grid_orog%nx, params%g_grid_orog%ny))
    allocate(alb_temp (params%g_grid%nx, params%g_grid%ny))
    allocate(if_temp  (params%g_grid%nx, params%g_grid%ny))
    allocate(vf_temp  (params%g_grid%nx, params%g_grid%ny))
    allocate(sif_temp (params%g_grid%nx, params%g_grid%ny))
    allocate(svf_temp (params%g_grid%nx, params%g_grid%ny))
    allocate(sd_temp  (params%g_grid%nx, params%g_grid%ny))

    if (GLC_DEBUG .and. main_task) then
       write(stdout,*) 'Upscale and splice the initial fields'
    end if

    ! Get initial fields from instances, splice together and return

    do i=1,params%ninstances

       call glint_upscaling(params%instances(i),   &
                            orog_temp,   alb_temp, &
                            if_temp,     vf_temp,  &
                            sif_temp,    svf_temp, &
                            sd_temp)

       ! Add this contribution to the global output
       ! Only the main task has valid values for the global output fields
       !
       ! TODO: Consider whether area_weighting should be true or false for these...
       ! arbitrarily setting them as true for now (I think that preserves their old
       ! behavior). But more generally: do we need the upscaling at all for this non-gcm
       ! case?
       if (main_task) then
          
          if (present(orog)) &
               orog = splice_field(orog, orog_temp, params%instances(i)%frac_cov_orog,&
                                       area_weighting=.true.)

          if (present(albedo)) &
               albedo = splice_field(albedo, alb_temp, params%instances(i)%frac_coverage,&
                                       area_weighting=.true.)

          if (present(ice_frac)) &
               ice_frac = splice_field(ice_frac, if_temp, params%instances(i)%frac_coverage,&
                                       area_weighting=.true.)

          if (present(veg_frac)) &
               veg_frac = splice_field(veg_frac, vf_temp, params%instances(i)%frac_coverage,&
                                       area_weighting=.true.)

          if (present(snowice_frac)) &
               snowice_frac = splice_field(snowice_frac,sif_temp,params%instances(i)%frac_coverage,&
                                       area_weighting=.true.)

          if (present(snowveg_frac)) &
               snowveg_frac = splice_field(snowveg_frac,svf_temp,params%instances(i)%frac_coverage,&
                                       area_weighting=.true.)

          if (present(snow_depth)) &
               snow_depth = splice_field(snow_depth,sd_temp,params%instances(i)%frac_coverage,&
                                       area_weighting=.true.)

       end if

    end do  ! ninstances

    ! Deallocate

    deallocate(orog_temp, alb_temp, if_temp, vf_temp, sif_temp, svf_temp,sd_temp)

    ! Sort out snow_model flag

    if (present(snow_model)) then
       snow_model = .false.
       do i=1, params%ninstances
          snow_model = (snow_model .or. glint_has_snow_model(params%instances(i)))
       end do
    end if

    ! Set output flag

    if (present(output_flag)) output_flag = .true.

  end subroutine initialise_glint

  !================================================================================

  subroutine initialise_glint_gcm(params,                         &
                                  lats,         longs,            &
                                  time_step,    paramfile,        &
                                  daysinyear,   start_time,       &
                                  ice_dt,       output_flag,      &
                                  glc_nec,                        &
                                  gfrac,        gtopo,            &
                                  grofi,        grofl,            &
                                  ice_sheet_grid_mask,            &
                                  icemask_coupled_fluxes,         &
                                  ghflx,        gmask,            &
                                  gcm_restart,  gcm_restart_file, &
                                  gcm_debug,    gcm_fileunit,     &
                                  test_coupling)

    ! Initialise the model for runs coupled to a GCM.
    ! For a multi-processor run, the main task should specify lats & longs spanning
    !  the full global domain; the other tasks should give 0-size lats & longs arrays
    ! Output arrays on the global grid are only valid on the main task
    !
    ! Note about ice_sheet_grid_mask and icemask_coupled_fluxes: ice_sheet_grid_mask is
    ! non-zero wherever CISM is operating - i.e., grid cells with icesheet or bare land
    ! (but not ocean). icemask_coupled_fluxes is similar, but is 0 for icesheet instances
    ! that have zero_gcm_fluxes = .true. Thus, icemask_coupled_fluxes can be used to
    ! determine the regions of the world in which CISM is operating and potentially
    ! sending non-zero fluxes to the climate model.

    use glimmer_config
    use glint_initialise
    use glimmer_log
    use glimmer_filenames
    use glimmer_physcon, only: rearth
    use glint_upscale, only: glint_upscaling_gcm
    use parallel, only: main_task

    implicit none

    ! Subroutine argument declarations --------------------------------------------------------

    type(glint_params),              intent(inout) :: params      !> parameters to be set
    real(dp),dimension(:),             intent(in)  :: lats,longs  !> location of gridpoints 
                                                                  !> in global data.
    integer,                           intent(in)  :: time_step   !> Timestep of calling model (hours)
    character(*),dimension(:),         intent(in)  :: paramfile   !> array of configuration filenames.
    integer,                  optional,intent(in)  :: daysinyear  !> Number of days in the year
    integer,                  optional,intent(in)  :: start_time  !> Time of first call to glint (hours)
    integer,                  optional,intent(out) :: ice_dt      !> Ice dynamics time-step in hours
    logical,                  optional,intent(out) :: output_flag !> Flag to show output set (provided for consistency)
    integer,                  optional,intent(in)  :: glc_nec     !> number of elevation classes for GCM input
    real(dp),dimension(:,:,0:),optional,intent(out) :: gfrac       !> ice+bare land fractional area [0,1]
    real(dp),dimension(:,:,0:),optional,intent(out) :: gtopo       !> surface elevation (m)
    real(dp),dimension(:,:,0:),optional,intent(out) :: ghflx       !> heat flux (W/m^2, positive down)
    real(dp),dimension(:,:),  optional,intent(out) :: grofi       !> ice runoff (kg/m^2/s = mm H2O/s)
    real(dp),dimension(:,:),  optional,intent(out) :: grofl       !> liquid runoff (kg/m^2/s = mm H2O/s)
    real(dp),dimension(:,:),  optional,intent(out) :: ice_sheet_grid_mask !mask of ice sheet grid coverage
    real(dp),dimension(:,:),  optional,intent(out) :: icemask_coupled_fluxes !mask of ice sheet grid coverage where we are potentially sending non-zero fluxes
    integer, dimension(:,:),  optional,intent(in)  :: gmask       !> mask = 1 where global data are valid
    logical,                  optional,intent(in)  :: gcm_restart ! logical flag to restart from a GCM restart file
    character(*),             optional,intent(in)  :: gcm_restart_file ! restart filename for a GCM restart
                                                                  ! (currently assumed to be CESM)
    logical,                  optional,intent(in)  :: gcm_debug   ! logical flag from GCM to output debug information
    integer,                  optional,intent(in)  :: gcm_fileunit! fileunit for reading config files
    logical,                  optional,intent(in)  :: test_coupling ! if true, force frequent coupling for testing purposes

    ! Internal variables -----------------------------------------------------------------------

    type(ConfigSection), pointer :: global_config, instance_config, section  ! configuration stuff

    character(len=100) :: message                                            ! For log-writing
    character(fname_length),dimension(:),pointer :: config_fnames=>null()    ! array of config filenames

    integer :: i, j

    integer,dimension(:),allocatable :: mbts,idts ! Array of mass-balance and ice dynamics timesteps

    real(dp),dimension(:,:,:),allocatable ::   &
               gfrac_temp, gtopo_temp, ghflx_temp ! Temporary output arrays

    real(dp),dimension(:,:),allocatable ::   &
               grofi_temp, grofl_temp, ice_sheet_grid_mask_temp, icemask_coupled_fluxes_temp  ! Temporary output arrays

    integer :: n
    integer :: nec       ! number of elevation classes

    if (present(gcm_debug)) then
       GLC_DEBUG = gcm_debug
    endif

    if (GLC_DEBUG .and. main_task) then
       write(stdout,*) 'Initializing glint'
    end if

    ! Initialise start time and calling model time-step (time_step = integer number of hours)
    ! We ignore t=0 by default 

    params%time_step = time_step

    ! Note: start_time = nhour_glint = 0 for an initial run.
    !       Does this create problems given that Glint convention is to ignore t = 0?

    if (present(start_time)) then
       params%start_time = start_time
    else
       params%start_time = time_step
    end if

    params%next_av_start = params%start_time

    params%gcm_restart = .false.
    if (present(gcm_restart)) then
       params%gcm_restart = gcm_restart
    endif

    params%gcm_restart_file = ''
    if (present(gcm_restart_file)) then
       params%gcm_restart_file = gcm_restart_file
    endif

    params%gcm_fileunit = 99
    if (present(gcm_fileunit)) then
       params%gcm_fileunit = gcm_fileunit
    endif

    nec = 1
    if (present(glc_nec)) then
       nec = glc_nec
    endif

    if (GLC_DEBUG .and. main_task) then
       write(stdout,*) 'time_step     =', params%time_step
       write(stdout,*) 'start_time    =', params%start_time
       write(stdout,*) 'next_av_start =', params%next_av_start
    end if

    ! Initialise year-length -------------------------------------------------------------------

    if (present(daysinyear)) then
       call glad_set_year_length(daysinyear)
    end if

    if (GLC_DEBUG .and. main_task) then
       write(stdout,*) 'Initialize global grid: present(gmask) =', present(gmask)
    end if

    ! Initialise main global grid --------------------------------------------------------------

    if (present(gmask)) then
       call new_global_grid(params%g_grid, longs, lats, nec=nec, mask=gmask, radius=rearth)
    else
       call new_global_grid(params%g_grid, longs, lats, nec=nec, radius=rearth)
    endif

    if (GLC_DEBUG .and. main_task) then
       write (stdout,*) ' ' 
       write (stdout,*) 'time_step (hr)  =', params%time_step
       write (stdout,*) 'start_time (hr) =', params%start_time
       write (stdout,*) 'Called new_global_grid '
       write (stdout,*) 'g_grid%nx =', params%g_grid%nx
       write (stdout,*) 'g_grid%ny =', params%g_grid%ny
       write (stdout,*) ' '
       write (stdout,*) 'g_grid%lons =', params%g_grid%lons
       write (stdout,*) ' '
       write (stdout,*) 'g_grid%lats =', params%g_grid%lats
       write (stdout,*) ' '
       write (stdout,*) 'g_grid%lon_bound =', params%g_grid%lon_bound
       write (stdout,*) ' '
       write (stdout,*) 'g_grid%lat_bound =', params%g_grid%lat_bound
       do j = 5, 10
          write (stdout,*)
          write (stdout,*) 'j, g_grid%mask =', j, params%g_grid%mask(:,j)
       enddo
    end if

    ! Allocate arrays -----------------------------------------------

    allocate(params%total_coverage(params%g_grid%nx, params%g_grid%ny))

    allocate(params%g_av_qsmb (params%g_grid%nx, params%g_grid%ny, 0:params%g_grid%nec))
    allocate(params%g_av_tsfc (params%g_grid%nx, params%g_grid%ny, 0:params%g_grid%nec))
    allocate(params%g_av_topo (params%g_grid%nx, params%g_grid%ny, 0:params%g_grid%nec))

    ! Initialise arrays ---------------------------------------------

    params%g_av_qsmb(:,:,:) = 0.d0
    params%g_av_tsfc(:,:,:) = 0.d0
    params%g_av_topo(:,:,:) = 0.d0

    ! ---------------------------------------------------------------
    ! Zero coverage maps and normalisation fields for main grid
    ! (Not using orography grid for GCM coupling)
    ! ---------------------------------------------------------------

    params%total_coverage = 0.d0

    if (GLC_DEBUG .and. main_task) then
       write(stdout,*) 'Read paramfile'
       write(stdout,*) 'paramfile =', paramfile
    end if

    ! ---------------------------------------------------------------
    ! Determine how many instances there are, according to what
    ! configuration files we've been provided with
    ! ---------------------------------------------------------------

    if (size(paramfile) == 1) then
       ! Load the configuration file into the linked list
       call ConfigRead(process_path(paramfile(1)), global_config, params%gcm_fileunit)    
       call glint_readconfig(global_config, params%ninstances, config_fnames, paramfile) ! Parse the list
    else
       params%ninstances = size(paramfile)
       allocate(config_fnames(params%ninstances))
       config_fnames = paramfile
    end if

    allocate(params%instances(params%ninstances))

    allocate(mbts(params%ninstances), idts(params%ninstances))

    if (GLC_DEBUG .and. main_task) then
       write(stdout,*) 'Number of instances =', params%ninstances
       write(stdout,*) 'Read config files and initialize each instance'
    end if

    ! ---------------------------------------------------------------
    ! Read config files, and initialise instances accordingly
    ! ---------------------------------------------------------------

    call write_log('Reading instance configurations')
    call write_log('-------------------------------')

    do i=1,params%ninstances

       call ConfigRead(process_path(config_fnames(i)),instance_config, params%gcm_fileunit)

         !WHL - I don't think this will be needed; commented out for now
!!       if (present(extraconfigs)) then
!!          if (size(extraconfigs)>=i) then
!!             call ConfigCombine(instance_config,extraconfigs(i))
!!          end if
!!       end if

       call glint_i_initialise_gcm(instance_config,     params%instances(i),     &
                                   params%g_grid,                                &
                                   mbts(i),             idts(i),                 &
                                   params%start_time,   params%time_step,        &
                                   params%gcm_restart,  params%gcm_restart_file, &
                                   params%gcm_fileunit, test_coupling )

       params%total_coverage = params%total_coverage + params%instances(i)%frac_coverage

    end do    ! ninstances

    ! Check that all mass-balance time-steps are the same length and
    ! assign that value to the top-level variable

    params%tstep_mbal = check_mbts(mbts)

    if (present(ice_dt)) then
       ice_dt = check_mbts(idts)
    end if

    if (GLC_DEBUG .and. main_task) then
       write(stdout,*) 'tstep_mbal =', params%tstep_mbal
       write(stdout,*) 'start_time =', params%start_time
       write(stdout,*) 'time_step =',  params%time_step
       if (present(ice_dt)) write(stdout,*) 'ice_dt =', ice_dt
    end if

    ! Check time-steps divide into one another appropriately.

    if (.not. (mod (params%tstep_mbal, params%time_step) == 0)) then
       call write_log('The mass-balance timestep must be an integer multiple of the forcing time-step', &
                       GM_FATAL,__FILE__,__LINE__)
    end if

    ! Make sure we don't have coverage greater than one at any point.

    where (params%total_coverage > 1.d0) params%total_coverage = 1.d0
    params%coverage_calculated = .true.

    ! Zero optional outputs, if present

    if (present(gfrac)) gfrac(:,:,:) = 0.d0
    if (present(gtopo)) gtopo(:,:,:) = 0.d0
    if (present(ghflx)) ghflx(:,:,:) = 0.d0
    if (present(grofi)) grofi(:,:)   = 0.d0
    if (present(grofl)) grofl(:,:)   = 0.d0
    if (present(ice_sheet_grid_mask))  ice_sheet_grid_mask(:,:)   = 0.d0
    if (present(icemask_coupled_fluxes))  icemask_coupled_fluxes(:,:)   = 0.d0

    ! Allocate arrays

    allocate(gfrac_temp(params%g_grid%nx, params%g_grid%ny, 0:params%g_grid%nec))
    allocate(gtopo_temp(params%g_grid%nx, params%g_grid%ny, 0:params%g_grid%nec))
    allocate(ghflx_temp(params%g_grid%nx, params%g_grid%ny, 0:params%g_grid%nec))
    allocate(grofi_temp(params%g_grid%nx, params%g_grid%ny))
    allocate(grofl_temp(params%g_grid%nx, params%g_grid%ny))
    allocate(ice_sheet_grid_mask_temp(params%g_grid%nx, params%g_grid%ny))
    allocate(icemask_coupled_fluxes_temp(params%g_grid%nx, params%g_grid%ny))

    if (GLC_DEBUG .and. main_task) then
       write(stdout,*) 'Upscale and splice the initial fields'
    end if

    ! Get initial fields from instances, splice together and return

    do i = 1, params%ninstances

       ! Upscale the output fields for this instance

       if (GLC_DEBUG .and. main_task) then
          print*, 'Do initial upscaling, i =', i
       endif

       call glint_upscaling_gcm(params%instances(i), params%g_grid%nec, &
                                params%instances(i)%lgrid%size%pt(1),   &
                                params%instances(i)%lgrid%size%pt(2),   &
                                params%g_grid%nx,   params%g_grid%ny,   &
                                params%g_grid%box_areas,                &   
                                gfrac_temp,          gtopo_temp,        &
                                grofi_temp,          grofl_temp,        &
                                ghflx_temp,                             &
                                init_call = .true.)

       call compute_ice_sheet_grid_mask(ice_sheet_grid_mask_temp, gfrac_temp)
       call compute_icemask_coupled_fluxes(icemask_coupled_fluxes_temp, &
                                           ice_sheet_grid_mask_temp, &
                                           params%instances(i))

       ! Splice together with the global output

       if (GLC_DEBUG .and. main_task) then
          print*, 'Spliced, i =', i
       endif

       call splice_fields_gcm(gfrac_temp, gtopo_temp,    &
                              grofi_temp, grofl_temp,    &
                              ghflx_temp,                &
                              ice_sheet_grid_mask_temp,  &
                              icemask_coupled_fluxes_temp, &
                              gfrac,      gtopo,         &
                              grofi,      grofl,         &
                              ghflx,                     &
                              ice_sheet_grid_mask,       &
                              icemask_coupled_fluxes,    &
                              params%g_grid%nec,         &
                              params%instances(i)%frac_coverage)

    end do       ! ninstances

    ! Deallocate

    deallocate(gfrac_temp, gtopo_temp, grofi_temp, grofl_temp, ghflx_temp)
    deallocate(ice_sheet_grid_mask_temp, icemask_coupled_fluxes_temp)

    ! Set output flag
    if (present(output_flag)) output_flag = .true.

    if (GLC_DEBUG .and. main_task) then
       write(stdout,*) 'Done in initialise_glint_gcm'
    endif

  end subroutine initialise_glint_gcm

  !================================================================================

  subroutine glint(params,         time,            &
                   rawtemp,        rawprecip,       &
                   orog,                            &
                   zonwind,        merwind,         &
                   humid,          lwdown,          &
                   swdown,         airpress,        &
                   output_flag,                     &
                   orog_out,       albedo,          &
                   ice_frac,       veg_frac,        &
                   snowice_frac,   snowveg_frac,    &
                   snow_depth,                      &
                   water_in,       water_out,       &
                   total_water_in, total_water_out, &
                   ice_volume,     ice_tstep)

    !> Main Glint subroutine.
    !>
    !> This should be called daily or hourly, depending on
    !> the mass-balance scheme being used. It does all necessary 
    !> spatial and temporal averaging, and calls the dynamics 
    !> part of the model when required. 
    !>
    !> Input fields should be taken as means over the period since the last call.
    !> See the user documentation for more information.
    !>
    !> Global output fields are only valid on the main task. Fields that are integrated
    !> over the whole domain (total_water_in, total_water_out, ice_volume) are only
    !> valid in single-task runs; trying to compute these in multi-task runs will generate a
    !> fatal error.
    !>
    !> Note that the total ice volume returned is the total at the end of the time-step;
    !> the water fluxes are valid over the duration of the timestep. Thus the difference
    !> between \texttt{total\_water\_in} and \texttt{total\_water\_out} should be equal
    !> to the change in \texttt{ice\_volume}, after conversion between m$^3$ and kg.

    use glimmer_utils
    use glint_interp
    use glint_timestep, only: glint_i_tstep
    use glint_downscale, only: glint_downscaling
    use glint_upscale, only: glint_upscaling
    use glimmer_log
    use glimmer_paramets, only: scyr
    use parallel, only: main_task, tasks
    implicit none

    ! Subroutine argument declarations -------------------------------------------------------------

    type(glint_params),              intent(inout) :: params          !> parameters for this run
    integer,                         intent(in)    :: time            !> Current model time        (hours)
    real(dp),dimension(:,:),target,  intent(in)    :: rawtemp         !> Surface temperature field (deg C)
    real(dp),dimension(:,:),target,  intent(in)    :: rawprecip       !> Precipitation rate        (mm/s)
    real(dp),dimension(:,:),         intent(in)    :: orog            !> The large-scale orography (m)
    real(dp),dimension(:,:),optional,intent(in)    :: zonwind,merwind !> Zonal and meridional components 
                                                                      !> of the wind field         (m/s)
    real(dp),dimension(:,:),optional,intent(in)    :: humid           !> Surface humidity (%)
    real(dp),dimension(:,:),optional,intent(in)    :: lwdown          !> Downwelling longwave (W/m$^2$)
    real(dp),dimension(:,:),optional,intent(in)    :: swdown          !> Downwelling shortwave (W/m$^2$)
    real(dp),dimension(:,:),optional,intent(in)    :: airpress        !> surface air pressure (Pa)
    logical,                optional,intent(out)   :: output_flag     !> Set true if outputs set
    real(dp),dimension(:,:),optional,intent(inout) :: orog_out        !> The fed-back, output orography (m)
    real(dp),dimension(:,:),optional,intent(inout) :: albedo          !> surface albedo
    real(dp),dimension(:,:),optional,intent(inout) :: ice_frac        !> grid-box ice-fraction
    real(dp),dimension(:,:),optional,intent(inout) :: veg_frac        !> grid-box veg-fraction
    real(dp),dimension(:,:),optional,intent(inout) :: snowice_frac    !> grid-box snow-covered ice fraction
    real(dp),dimension(:,:),optional,intent(inout) :: snowveg_frac    !> grid-box snow-covered veg fraction
    real(dp),dimension(:,:),optional,intent(inout) :: snow_depth      !> grid-box mean snow depth (m water equivalent)
    real(dp),dimension(:,:),optional,intent(inout) :: water_in        !> Input water flux          (mm)
    real(dp),dimension(:,:),optional,intent(inout) :: water_out       !> Output water flux         (mm)
    real(dp),               optional,intent(inout) :: total_water_in  !> Area-integrated water flux in (kg)
    real(dp),               optional,intent(inout) :: total_water_out !> Area-integrated water flux out (kg)
    real(dp),               optional,intent(inout) :: ice_volume      !> Total ice volume (m$^3$)
    logical,                optional,intent(out)   :: ice_tstep       !> Set when an ice-timestep has been done, and
                                                                      !> water balance information is available

    ! Internal variables ----------------------------------------------------------------------------

    integer :: i, j
    real(dp),dimension(:,:),allocatable :: albedo_temp, if_temp, vf_temp, sif_temp, svf_temp,  &
                                           sd_temp, wout_temp, orog_out_temp, win_temp
    real(dp) :: twin_temp,twout_temp,icevol_temp
    type(output_flags) :: out_f
    logical :: icets
    character(250) :: message
    real(dp),dimension(size(rawprecip,1),size(rawprecip,2)),target :: anomprecip
    real(dp),dimension(size(rawtemp,1),  size(rawtemp,2)),  target :: anomtemp
    real(dp),dimension(:,:),pointer :: precip
    real(dp),dimension(:,:),pointer :: temp
    real(dp) :: yearfrac

    if (GLC_DEBUG .and. main_task) then
!       write (stdout,*) 'In subroutine glint, current time (hr) =', time
!       write (stdout,*) 'av_start_time =', params%av_start_time
!       write (stdout,*) 'next_av_start =', params%next_av_start
!       write (stdout,*) 'new_av =', params%new_av
!       write (stdout,*) 'tstep_mbal =', params%tstep_mbal
    end if

    ! Check we're expecting a call now --------------------------------------------------------------

    if (params%new_av) then
       if (time == params%next_av_start) then
          params%av_start_time = time
          params%new_av = .false.
       else
          write(message,*) 'Unexpected calling of GLINT at time ', time
          call write_log(message,GM_FATAL,__FILE__,__LINE__)
       end if
    else
       if (mod(time-params%av_start_time,params%time_step) /= 0) then
          write(message,*) 'Unexpected calling of GLINT at time ', time
          call write_log(message,GM_FATAL,__FILE__,__LINE__)
       end if
    end if

    ! Check input fields are correct ----------------------------------------------------------------

    call check_input_fields(params, humid, lwdown, swdown, airpress, zonwind, merwind)

    ! Reset output flag

    if (present(output_flag)) output_flag = .false.
    if (present(ice_tstep))   ice_tstep = .false.

    ! Sort out anomaly coupling

    if (params%anomaly_params%enabled) then
       yearfrac = real(mod(time,days_in_year),dp)/real(days_in_year,dp)
       call anomaly_calc(params%anomaly_params, yearfrac, rawtemp, rawprecip, anomtemp, anomprecip)
       precip => anomprecip
       temp   => anomtemp
    else
       precip => rawprecip
       temp   => rawtemp
    end if

    ! Do averaging and so on...

    call accumulate_averages(params,           &
                             temp,    precip,  &
                             zonwind, merwind, &
                             humid,   lwdown,  &
                             swdown,  airpress)

    ! Increment step counter

    params%av_steps = params%av_steps + 1

    ! ---------------------------------------------------------
    ! If this is a mass balance timestep, prepare global fields, and do a timestep
    ! for each model instance
    ! ---------------------------------------------------------

    if (time - params%av_start_time + params%time_step > params%tstep_mbal) then

       write(message,*) &
            'Incomplete forcing of GLINT mass-balance time-step detected at time ', time
       call write_log(message,GM_FATAL,__FILE__,__LINE__)

    else if (time - params%av_start_time + params%time_step == params%tstep_mbal) then

       ! Set output_flag

       ! At present, outputs are done for each mass-balance timestep, since
       ! that involved least change to the code. However, it might be good
       ! to change the output to occur with user-specified frequency.

       if (present(output_flag)) output_flag = .true.

       ! Allocate output fields

       if (present(orog_out)) then
          allocate(orog_out_temp(size(orog_out,1),size(orog_out,2)))
       else
          allocate(orog_out_temp(params%g_grid_orog%nx, params%g_grid_orog%ny))
       end if
       allocate(albedo_temp(size(orog,1),size(orog,2)))
       allocate(if_temp(size(orog,1),size(orog,2)))
       allocate(vf_temp(size(orog,1),size(orog,2)))
       allocate(sif_temp(size(orog,1),size(orog,2)))
       allocate(svf_temp(size(orog,1),size(orog,2)))
       allocate(sd_temp(size(orog,1),size(orog,2)))
       allocate(wout_temp(size(orog,1),size(orog,2)))
       allocate(win_temp(size(orog,1),size(orog,2)))

       ! Populate output flag derived type

       call populate_output_flags(out_f,                           &
                                  orog_out,       albedo,          &
                                  ice_frac,       veg_frac,        &
                                  snowice_frac,   snowveg_frac,    &
                                  snow_depth,                      &
                                  water_in,       water_out,       &
                                  total_water_in, total_water_out, &
                                  ice_volume)

       ! Zero outputs if present

       if (present(orog_out))        orog_out        = 0.d0
       if (present(albedo))          albedo          = 0.d0
       if (present(ice_frac))        ice_frac        = 0.d0
       if (present(veg_frac))        veg_frac        = 0.d0
       if (present(snowice_frac))    snowice_frac    = 0.d0
       if (present(snowveg_frac))    snowveg_frac    = 0.d0
       if (present(snow_depth))      snow_depth      = 0.d0
       if (present(water_out))       water_out       = 0.d0
       if (present(water_in))        water_in        = 0.d0
       if (present(total_water_in))  total_water_in  = 0.d0
       if (present(total_water_out)) total_water_out = 0.d0
       if (present(ice_volume))      ice_volume      = 0.d0

       ! Calculate averages by dividing by number of steps elapsed
       ! since last model timestep.

       call calculate_averages(params)

       ! Calculate total accumulated precipitation - multiply
       ! by time since last model timestep

       params%g_av_precip = params%g_av_precip*params%tstep_mbal*hours2seconds

       ! Calculate temperature half-range

       params%g_temp_range = (params%g_max_temp-params%g_min_temp)/2.0

       write(stdout,*) 'Take a mass balance timestep, time (hr) =', time
       write(stdout,*) 'av_steps =', real(params%av_steps,dp)
       write(stdout,*) 'tstep_mbal (hr) =', params%tstep_mbal

       ! Do a timestep for each instance

       do i = 1,params%ninstances
          
          !WHL - Moved some code here from start of glint_i_tstep

          ! Check whether we're doing anything this time   

          if (time == params%instances(i)%next_time) then

             params%instances(i)%next_time = params%instances(i)%next_time + params%instances(i)%mbal_tstep

             ! Downscale input fields from global to local grid

             call glint_downscaling(params%instances(i),                           &
                                    params%g_av_temp,     params%g_temp_range,     &
                                    params%g_av_precip,   orog,                    &
                                    params%g_av_zonwind,  params%g_av_merwind,     &
                                    params%g_av_humid,    params%g_av_lwdown,      &
                                    params%g_av_swdown,   params%g_av_airpress,    &
                                    .true.)   ! orogflag = .true.

             call glint_i_tstep(time,                    params%instances(i),       &
                                orog_out_temp,                                      &
                                albedo_temp,             if_temp,                   &
                                vf_temp,                 sif_temp,                  &
                                svf_temp,                sd_temp,                   &
                                win_temp,                wout_temp,                 &
                                twin_temp,               twout_temp,                &
                                icevol_temp,             out_f,                     &
                                icets)

             if (GLC_DEBUG .and. main_task) then
                write(stdout,*) 'Finished glc_glint_ice tstep, instance =', i
                write(stdout,*) 'Upscale fields to global grid'
             end if

             ! Add this contribution to the global output
             ! Only the main task has valid values for the global output fields
             !
             ! TODO: Consider whether area_weighting should be true or false for these...
             ! arbitrarily setting them as true for now (I think that preserves their old
             ! behavior). But more generally: do we need the upscaling at all for this
             ! non-gcm case?
             if (main_task) then
             
                if (present(orog_out)) &
                     orog_out = splice_field(orog_out,orog_out_temp, &
                                       params%instances(i)%frac_cov_orog,&
                                       area_weighting=.true.)

                if (present(albedo)) &
                     albedo = splice_field(albedo,albedo_temp, &
                                       params%instances(i)%frac_coverage,&
                                       area_weighting=.true.)

                if (present(ice_frac)) &
                     ice_frac = splice_field(ice_frac,if_temp, &
                                       params%instances(i)%frac_coverage,&
                                       area_weighting=.true.)

                if (present(veg_frac)) &
                     veg_frac = splice_field(veg_frac,vf_temp, &
                                       params%instances(i)%frac_coverage,&
                                       area_weighting=.true.)

                if (present(snowice_frac)) &
                     snowice_frac = splice_field(snowice_frac,sif_temp, &
                                       params%instances(i)%frac_coverage,&
                                       area_weighting=.true.)

                if (present(snowveg_frac)) &
                     snowveg_frac = splice_field(snowveg_frac, &
                                       svf_temp,params%instances(i)%frac_coverage,&
                                       area_weighting=.true.)

                if (present(snow_depth)) &
                     snow_depth = splice_field(snow_depth, &
                                       sd_temp,params%instances(i)%frac_coverage,&
                                       area_weighting=.true.)

                if (present(water_in)) &
                     water_in = splice_field(water_in,win_temp, &
                                       params%instances(i)%frac_coverage,&
                                       area_weighting=.true.)

                if (present(water_out)) &
                     water_out = splice_field(water_out, &
                                       wout_temp, params%instances(i)%frac_coverage,&
                                       area_weighting=.true.)

             end if

             ! Add total water variables to running totals
             ! WJS (1-15-13): These fields are only valid in single-task runs; multi-task
             ! runs should generate an error in glint_i_tstep if you try to compute any of
             ! these. But to be safe, we check here, too

             if (present(total_water_in))  then
                if (tasks > 1) call write_log('total_water_in is only valid when running with a single task', &
                                              GM_FATAL, __FILE__, __LINE__)

                total_water_in  = total_water_in  + twin_temp
             end if

             if (present(total_water_out)) then
                if (tasks > 1) call write_log('total_water_out is only valid when running with a single task', &
                                              GM_FATAL, __FILE__, __LINE__)

                total_water_out = total_water_out + twout_temp
             end if

             if (present(ice_volume)) then
                if (tasks > 1) call write_log('ice_volume is only valid when running with a single task', &
                                              GM_FATAL, __FILE__, __LINE__)

                ice_volume      = ice_volume      + icevol_temp
             end if

             ! Set flag
             if (present(ice_tstep)) then
                ice_tstep = (ice_tstep .or. icets)
             end if

          endif ! time = next_time

       enddo    ! ninstances

       ! Scale output water fluxes to be in mm/s

       if (present(water_in))  water_in  = water_in/ &
                                           (params%tstep_mbal*hours2seconds)

       if (present(water_out)) water_out = water_out/ &
                                           (params%tstep_mbal*hours2seconds)

       ! ---------------------------------------------------------
       ! Reset averaging fields, flags and counters
       ! ---------------------------------------------------------

       params%g_av_temp    = 0.d0
       params%g_av_precip  = 0.d0
       params%g_av_zonwind = 0.d0
       params%g_av_merwind = 0.d0
       params%g_av_humid   = 0.d0
       params%g_av_lwdown  = 0.d0
       params%g_av_swdown  = 0.d0
       params%g_av_airpress= 0.d0
       params%g_temp_range = 0.d0
       params%g_max_temp   = -1000.d0
       params%g_min_temp   =  1000.d0

       params%av_steps      = 0
       params%new_av        = .true.
       params%next_av_start = time+params%time_step

       deallocate(albedo_temp,if_temp,vf_temp,sif_temp,svf_temp,sd_temp,wout_temp,win_temp,orog_out_temp)

    endif    ! time - params%av_start_time + params%time_step = params%tstep_mbal

  end subroutine glint

  !===================================================================

  subroutine glint_gcm(params,         time,            &
                       qsmb,           tsfc,            &
                       topo,                            &
                       output_flag,    ice_tstep,       &
                       gfrac,          gtopo,           &
                       grofi,          grofl,           &
                       ice_sheet_grid_mask,             &
                       icemask_coupled_fluxes,          &
                       ghflx)

    ! Main Glint subroutine for GCM coupling.
    !
    ! It does all necessary spatial and temporal averaging, 
    ! and calls the dynamic ice sheet model when required. 
    !
    ! Input fields should be taken as means over the period since the last call.
    ! See the user documentation for more information.
    !
    ! Global output fields are only valid on the main task. Fields that are integrated
    ! over the whole domain are only valid in single-task runs; trying to compute these 
    ! in multi-task runs will generate a fatal error.
    !
    ! Note about ice_sheet_grid_mask and icemask_coupled_fluxes: ice_sheet_grid_mask is
    ! non-zero wherever CISM is operating - i.e., grid cells with icesheet or bare land
    ! (but not ocean). icemask_coupled_fluxes is similar, but is 0 for icesheet instances
    ! that have zero_gcm_fluxes = .true. Thus, icemask_coupled_fluxes can be used to
    ! determine the regions of the world in which CISM is operating and potentially
    ! sending non-zero fluxes to the climate model.

    use glimmer_utils
    use glint_interp
    use glint_timestep, only: glint_i_tstep_gcm
    use glint_downscale, only: glint_downscaling_gcm
    use glint_upscale, only: glint_upscaling_gcm
    use glimmer_log
    use glimmer_paramets, only: scyr
    use parallel, only: main_task, tasks

    implicit none

    ! Subroutine argument declarations -------------------------------------------------------------

    type(glint_params),              intent(inout) :: params          !> parameters for this run
    integer,                         intent(in)    :: time            !> Current model time        (hours)

    real(dp),dimension(:,:,0:),intent(in)    :: qsmb          ! input surface mass balance of glacier ice (kg/m^2/s)
    real(dp),dimension(:,:,0:),intent(in)    :: tsfc          ! input surface ground temperature (deg C)
    real(dp),dimension(:,:,0:),intent(in)    :: topo          ! input surface elevation (m)

    logical,                optional,intent(out)   :: output_flag     ! Set true if outputs are set
    logical,                optional,intent(out)   :: ice_tstep       ! Set when an ice dynamic timestep has been done
                                                                      !  and new output is available

    real(dp),dimension(:,:,:),optional,intent(inout) :: gfrac         ! output ice fractional area [0,1]
    real(dp),dimension(:,:,:),optional,intent(inout) :: gtopo         ! output surface elevation (m)
    real(dp),dimension(:,:,:),optional,intent(inout) :: ghflx         ! output heat flux (W/m^2, positive down)
    real(dp),dimension(:,:),  optional,intent(inout) :: grofi         ! output ice runoff (kg/m^2/s = mm H2O/s)
    real(dp),dimension(:,:),  optional,intent(inout) :: grofl         ! output liquid runoff (kg/m^2/s = mm H2O/s)
    real(dp),dimension(:,:),  optional,intent(inout) :: ice_sheet_grid_mask !mask of ice sheet grid coverage
    real(dp),dimension(:,:),  optional,intent(inout) :: icemask_coupled_fluxes !mask of ice sheet grid coverage where we are potentially sending non-zero fluxes

    ! Internal variables ----------------------------------------------------------------------------

    integer :: i

    logical :: icets
    character(250) :: message

    real(dp), dimension(:,:,:), allocatable ::   &
       gfrac_temp    ,&! gfrac for a single instance
       gtopo_temp    ,&! gtopo for a single instance
       ghflx_temp      ! ghflx for a single instance

    real(dp), dimension(:,:), allocatable ::   &
       grofi_temp    ,&! grofi for a single instance
       grofl_temp    ,&! grofl for a single instance
       ice_sheet_grid_mask_temp, &    ! ice_sheet_grid_mask for a single instance
       icemask_coupled_fluxes_temp    ! icemask_coupled_fluxes for a single instance

    if (GLC_DEBUG .and. main_task) then
       if (params%new_av) then
          write (stdout,*) 'In subroutine glint_gcm, current time (hr) =', time
          write (stdout,*) 'av_start_time =', params%av_start_time
          write (stdout,*) 'next_av_start =', params%next_av_start
          write (stdout,*) 'new_av =', params%new_av
          write (stdout,*) 'tstep_mbal =', params%tstep_mbal
       endif
    end if

    ! Check we're expecting a call now --------------------------------------------------------------

    if (params%new_av) then
       if (time == params%next_av_start) then
          params%av_start_time = time
          params%new_av = .false.
       else
          write(message,*) 'Unexpected calling of GLINT at time ', time
          call write_log(message,GM_FATAL,__FILE__,__LINE__)
       end if
    else
       if (mod (time - params%av_start_time, params%time_step) /= 0) then
          write(message,*) 'Unexpected calling of GLINT at time ', time
          call write_log(message,GM_FATAL,__FILE__,__LINE__)
       end if
    end if

    ! Check input fields are correct ----------------------------------------------------------------

    ! Reset output flag

    if (present(output_flag)) output_flag = .false.
    if (present(ice_tstep))   ice_tstep = .false.

    ! Accumulate input fields for later averaging

    call accumulate_averages_gcm(params, qsmb, tsfc, topo)

    ! Increment step counter

    params%av_steps = params%av_steps + 1

    ! ---------------------------------------------------------
    ! If this is a mass balance timestep, prepare global fields, and do a timestep
    ! for each model instance
    ! ---------------------------------------------------------

    if (time - params%av_start_time + params%time_step > params%tstep_mbal) then

       write(message,*) &
            'Incomplete forcing of GLINT mass-balance time-step detected at time ', time
       call write_log(message,GM_FATAL,__FILE__,__LINE__)

    else if (time - params%av_start_time + params%time_step == params%tstep_mbal) then

       ! Set output_flag

       ! At present, outputs are done for each mass-balance timestep, since
       ! that involved least change to the code. However, it might be good
       ! to change the output to occur with user-specified frequency.

       if (present(output_flag)) output_flag = .true.

       ! Allocate output fields
       ! Each *_temp field contains the output for one ice sheet instance.
       ! If there are multiple instances, the various *_temp fields are spliced together.

       allocate(gfrac_temp(params%g_grid%nx, params%g_grid%ny, 0:params%g_grid%nec))
       allocate(gtopo_temp(params%g_grid%nx, params%g_grid%ny, 0:params%g_grid%nec))
       allocate(ghflx_temp(params%g_grid%nx, params%g_grid%ny, 0:params%g_grid%nec))
       allocate(grofi_temp(params%g_grid%nx, params%g_grid%ny))
       allocate(grofl_temp(params%g_grid%nx, params%g_grid%ny))
       allocate(ice_sheet_grid_mask_temp(params%g_grid%nx, params%g_grid%ny))
       allocate(icemask_coupled_fluxes_temp(params%g_grid%nx, params%g_grid%ny))

       ! Zero global outputs if present

       if (present(gfrac)) gfrac(:,:,:) = 0.d0
       if (present(gtopo)) gtopo(:,:,:) = 0.d0
       if (present(ghflx)) ghflx(:,:,:) = 0.d0
       if (present(grofi)) grofi(:,:)   = 0.d0
       if (present(grofl)) grofl(:,:)   = 0.d0
       if (present(ice_sheet_grid_mask))  ice_sheet_grid_mask(:,:)   = 0.d0
       if (present(icemask_coupled_fluxes))  icemask_coupled_fluxes(:,:)   = 0.d0

       ! Calculate averages by dividing by number of steps elapsed
       ! since last model timestep.

       call calculate_averages_gcm(params)

       ! Calculate total surface mass balance - multiply by time since last model timestep
       ! Note on units: We want g_av_qsmb to have units of meters w.e. (accumulated over mass balance time step)
       ! Initial units are kg m-2 s-1 = mm s-1
       ! Divide by 1000 to convert from mm to m
       ! Multiply by hours2seconds = 3600 to convert from 1/s to 1/hr.  (tstep_mbal has units of hours)

       !TODO - Modify code so that qsmb and acab are always in kg m-2 s-1 water equivalent?
       params%g_av_qsmb(:,:,:) = params%g_av_qsmb(:,:,:) * params%tstep_mbal * hours2seconds / 1000.d0

       ! Do a timestep for each instance

       do i = 1, params%ninstances
          
          if (time == params%instances(i)%next_time) then

             params%instances(i)%next_time = params%instances(i)%next_time + params%instances(i)%mbal_tstep

             ! Downscale input fields from global to local grid
             ! This subroutine computes instance%acab and instance%artm, the key inputs to Glide.

             if (GLC_DEBUG .and. main_task) write(stdout,*) 'Downscale fields to local grid, time (hr) =', time
             call glint_downscaling_gcm (params%instances(i),   &
                                         params%g_av_qsmb,      &
                                         params%g_av_tsfc,      &
                                         params%g_av_topo,      &
                                         params%g_grid%mask)

             if (GLC_DEBUG .and. main_task) write(stdout,*) 'Take a glint time step, instance', i
             call glint_i_tstep_gcm(time,                  &
                                    params%instances(i),   &
                                    icets)


             ! Set flag
             if (present(ice_tstep)) then
                ice_tstep = (ice_tstep .or. icets)
             end if

             ! Upscale the output to elevation classes on the global grid

             if (GLC_DEBUG .and. main_task) write(stdout,*) 'Upscale fields to global grid, time(hr) =', time
             call glint_upscaling_gcm(params%instances(i), params%g_grid%nec, &
                                      params%instances(i)%lgrid%size%pt(1),   &
                                      params%instances(i)%lgrid%size%pt(2),   &
                                      params%g_grid%nx,   params%g_grid%ny,   &                               
                                      params%g_grid%box_areas,                &
                                      gfrac_temp,          gtopo_temp,        &
                                      grofi_temp,          grofl_temp,        &
                                      ghflx_temp )
             
             call compute_ice_sheet_grid_mask(ice_sheet_grid_mask_temp, gfrac_temp)
             call compute_icemask_coupled_fluxes(icemask_coupled_fluxes_temp, &
                                                 ice_sheet_grid_mask_temp, &
                                                 params%instances(i))

             ! Add the contribution from this instance to the global output

             call splice_fields_gcm(gfrac_temp, gtopo_temp,    & !gfrac_temp here is fractional area, for each elevation level,
                                    grofi_temp, grofl_temp,    & !  of the total land+ice area
                                    ghflx_temp,                &
                                    ice_sheet_grid_mask_temp,  &
                                    icemask_coupled_fluxes_temp, &
                                    gfrac,      gtopo,         & !gfrac here is the fractional area, for each elevation level, 
                                    grofi,      grofl,         & ! of the fractional area of the total grid cell that is covered by CISM-owned land
                                    ghflx,                     &
                                    ice_sheet_grid_mask,       &
                                    icemask_coupled_fluxes,    &
                                    params%g_grid%nec,         &
                                    params%instances(i)%frac_coverage)
                         
          endif   ! time = next_time

       enddo    ! ninstances

       ! ---------------------------------------------------------
       ! Reset averaging fields, flags and counters
       ! ---------------------------------------------------------

       params%g_av_qsmb(:,:,:) = 0.d0
       params%g_av_tsfc(:,:,:) = 0.d0
       params%g_av_topo(:,:,:) = 0.d0

       params%av_steps      = 0
       params%new_av        = .true.
       params%next_av_start = time + params%time_step

       deallocate(gfrac_temp, gtopo_temp, grofi_temp, grofl_temp, ghflx_temp)
       deallocate(ice_sheet_grid_mask_temp, icemask_coupled_fluxes_temp)

       if (GLC_DEBUG .and. main_task) then
          write(stdout,*) 'Done in glint_gcm'
       endif

   endif    ! time - params%av_start_time + params%time_step > params%tstep_mbal

  end subroutine glint_gcm

  !===================================================================

  subroutine end_glint(params,close_logfile)

    !> tidy-up operations for Glint
    use glint_initialise
    use glimmer_log
    implicit none

    type(glint_params),intent(inout) :: params          ! parameters for this run
    logical, intent(in), optional    :: close_logfile   ! if true, then close the log file
                                                        ! (GCM may do this elsewhere)                                  
    integer :: i

    ! end individual instances

    do i = 1, params%ninstances
       call glint_i_end(params%instances(i))
    enddo

    if (present(close_logfile)) then
       if (close_logfile) call close_log
    else
       call close_log
    endif

  end subroutine end_glint

  !=====================================================

  integer function glint_coverage_map(params, coverage, cov_orog)

    !> Retrieve ice model fractional coverage map. 
    !> This function is provided so that glimmer may
    !> be restructured without altering the interface.
    !> This is currently only valid on the main task.
    !*RV Three return values are possible:
    !*RV \begin{description}
    !*RV \item[0 ---] Successful return
    !*RV \item[1 ---] Coverage map not calculated yet (fail)
    !*RV \item[2 ---] Coverage array is the wrong size (fail)
    !*RV \end{description}

    implicit none

    type(glint_params),intent(in) :: params                  !> ice model parameters
    real(dp),dimension(:,:),intent(out) :: coverage          !> array to hold coverage map
    real(dp),dimension(:,:),intent(out),optional :: cov_orog !> Orography coverage

    if (.not. params%coverage_calculated) then
       glint_coverage_map = 1
       return
    endif

    if ( size(coverage,1) /= params%g_grid%nx .or. &
         size(coverage,2) /= params%g_grid%ny) then
       glint_coverage_map = 2
       return
    endif

    glint_coverage_map = 0
    coverage = params%total_coverage
    if (present(cov_orog)) cov_orog = params%total_cov_orog

  end function glint_coverage_map

  !=====================================================

  !----------------------------------------------------------------------
  ! PRIVATE INTERNAL GLIMMER SUBROUTINES FOLLOW.............
  !----------------------------------------------------------------------

  subroutine glint_allocate_arrays(params)

    !> allocates glimmer arrays

    implicit none

    type(glint_params),intent(inout) :: params !> ice model parameters

    allocate(params%g_av_precip (params%g_grid%nx, params%g_grid%ny))
    allocate(params%g_av_temp   (params%g_grid%nx, params%g_grid%ny))
    allocate(params%g_max_temp  (params%g_grid%nx, params%g_grid%ny))
    allocate(params%g_min_temp  (params%g_grid%nx, params%g_grid%ny))
    allocate(params%g_temp_range(params%g_grid%nx, params%g_grid%ny))
    allocate(params%g_av_zonwind(params%g_grid%nx, params%g_grid%ny))
    allocate(params%g_av_merwind(params%g_grid%nx, params%g_grid%ny))
    allocate(params%g_av_humid  (params%g_grid%nx, params%g_grid%ny))
    allocate(params%g_av_lwdown (params%g_grid%nx, params%g_grid%ny))
    allocate(params%g_av_swdown (params%g_grid%nx, params%g_grid%ny))
    allocate(params%g_av_airpress(params%g_grid%nx,params%g_grid%ny))

    allocate(params%total_coverage(params%g_grid%nx, params%g_grid%ny))
    allocate(params%total_cov_orog(params%g_grid_orog%nx, params%g_grid_orog%ny))

  end subroutine glint_allocate_arrays

  !========================================================

  subroutine splice_fields_gcm(gfrac_temp, gtopo_temp,    &
                               grofi_temp, grofl_temp,    &
                               ghflx_temp,                &
                               ice_sheet_grid_mask_temp,  &
                               icemask_coupled_fluxes_temp, &
                               gfrac,      gtopo,         &
                               grofi,      grofl,         &
                               ghflx,                     &
                               ice_sheet_grid_mask,       &
                               icemask_coupled_fluxes,    &
                               nec,                       &
                               frac_coverage)

     use parallel, only: main_task

     ! Add the output for this instance to the global output

     real(dp), dimension(:,:,0:), intent(in) :: gfrac_temp  ! output fields for this instance
     real(dp), dimension(:,:,0:), intent(in) :: gtopo_temp  ! output fields for this instance
     real(dp), dimension(:,:,0:), intent(in) :: ghflx_temp  ! output fields for this instance
     real(dp), dimension(:,:),   intent(in) :: grofi_temp  ! output fields for this instance
     real(dp), dimension(:,:),   intent(in) :: grofl_temp  ! output fields for this instance
     real(dp), dimension(:,:),   intent(in) :: ice_sheet_grid_mask_temp  ! output fields for this instance
     real(dp), dimension(:,:),   intent(in) :: icemask_coupled_fluxes_temp  ! output fields for this instance

     real(dp), dimension(:,:,0:), intent(inout) :: gfrac    ! spliced global output field
     real(dp), dimension(:,:,0:), intent(inout) :: gtopo    ! spliced global output field
     real(dp), dimension(:,:,0:), intent(inout) :: ghflx    ! spliced global output field
     real(dp), dimension(:,:),   intent(inout) :: grofi    ! spliced global output field
     real(dp), dimension(:,:),   intent(inout) :: grofl    ! spliced global output field
     real(dp), dimension(:,:),   intent(inout) :: ice_sheet_grid_mask    ! spliced global output field
     real(dp), dimension(:,:),   intent(inout) :: icemask_coupled_fluxes ! spliced global output field

     integer, intent(in) :: nec   ! number of elevation classes

     real(dp), dimension(:,:), intent(in) :: frac_coverage  ! map of fractional coverage of global gridcells
                                                            ! by local gridcells
     
     integer :: n

     ! Only the main task has valid values for the global output fields

     if (main_task) then

        do n = 0, nec

           gfrac(:,:,n) = splice_field(gfrac(:,:,n),            &
                                       gfrac_temp(:,:,n),       &
                                       frac_coverage,           &
                                       area_weighting=.true.)

           gtopo(:,:,n) = splice_field(gtopo(:,:,n),            &
                                       gtopo_temp(:,:,n),       &
                                       frac_coverage,           &
                                       area_weighting=.false.)

           ! TODO: No thought has been given to whether area_weighting should be true or false for ghflx... 
           !       This needs to be considered once ghflx is hooked up to CLM.

           ghflx(:,:,n) = splice_field(ghflx(:,:,n),            &
                                       ghflx_temp(:,:,n),       &
                                       frac_coverage,           &
                                       area_weighting=.true.)

        enddo   ! nec

        ! area_weighting is false for grofi and grofl, because they are computed as sums
        ! per unit area of the GLOBAL grid (as opposed to averages per unit area of the
        ! icesheet grid). So the normalization done by the area_weighting option has, in
        ! effect, already been done.

        grofi(:,:) = splice_field(grofi(:,:),              &
                                  grofi_temp(:,:),         &
                                  frac_coverage,           &
                                  area_weighting=.false.)

        grofl(:,:) = splice_field(grofl(:,:),              &
                                  grofl_temp(:,:),         &
                                  frac_coverage,           &
                                  area_weighting=.false.)

        ! area_weighting for ice_sheet_grid_mask agrees with area_weighting for gfrac,
        ! since they are similar variables

        ice_sheet_grid_mask(:,:) = splice_field(ice_sheet_grid_mask(:,:), &
                                                ice_sheet_grid_mask_temp(:,:), &
                                                frac_coverage, &
                                                area_weighting = .true.)

        icemask_coupled_fluxes(:,:) = splice_field(icemask_coupled_fluxes(:,:), &
                                                   icemask_coupled_fluxes_temp(:,:), &
                                                   frac_coverage, &
                                                   area_weighting = .true.)

     endif  ! main_task

  end subroutine splice_fields_gcm

  !========================================================

  function splice_field(global, local, coverage, area_weighting)

    !> Splices an upscaled field into a global field

    ! Note that this does not handle multiple overlapping ice sheet instances

    real(dp),dimension(:,:),intent(in) :: global    !> Field to receive the splice
    real(dp),dimension(:,:),intent(in) :: local     !> The field to be spliced in
    real(dp),dimension(:,:),intent(in) :: coverage  !> The coverage fraction of the ice sheet grid on the global grid cell
    logical,intent(in) :: area_weighting            !> Do/not do area weighting
    real(dp),dimension(size(global,1),size(global,2)) :: splice_field

    where (coverage == 0.d0) splice_field = global
    
    if (area_weighting) then 
       where (coverage > 0.d0) 
          !Spliced field = area-weighted average of the local field and existing global field.
          splice_field = (global*(1.d0-coverage)) + (local*coverage)
       end where
    else
       where (coverage > 0.d0)
          !Spliced field is straight addition of local to global
          splice_field = global+local
       end where
    endif

  end function splice_field

  !========================================================

  !TODO - Move subroutine glint_readconfig to a glint_setup module, in analogy to glide_setup?

  subroutine glint_readconfig(config, ninstances, fnames, infnames)

    !> Determine whether a given config file is a
    !> top-level glint config file, and return parameters
    !> accordingly.

    use glimmer_config
    use glimmer_log
    implicit none

    ! Arguments -------------------------------------------

    type(ConfigSection),      pointer :: config !> structure holding sections of configuration file
    integer,              intent(out) :: ninstances !> Number of instances to create
    character(fname_length),dimension(:),pointer :: fnames !> list of filenames (output)
    character(fname_length),dimension(:) :: infnames !> list of filenames (input)

    ! Internal variables ----------------------------------

    type(ConfigSection), pointer :: section
    character(len=100) :: message
    integer :: i

    if (associated(fnames)) nullify(fnames)

    call GetSection(config,section,'GLINT')
    if (associated(section)) then
       call GetValue(section,'n_instance',ninstances)
       allocate(fnames(ninstances))
       do i=1,ninstances
          call GetSection(section%next,section,'GLINT instance')
          if (.not.associated(section)) then
             write(message,*) 'Must specify ',ninstances,' instance config files'
             call write_log(message,GM_FATAL,__FILE__,__LINE__)
          end if
          call GetValue(section,'name',fnames(i))
       end do
    else
       ninstances=1
       allocate(fnames(1))
       fnames=infnames
    end if

    ! Print some configuration information

!!$    call write_log('GLINT global')
!!$    call write_log('------------')
!!$    write(message,*) 'number of instances :',params%ninstances
!!$    call write_log(message)
!!$    call write_log('')

  end subroutine glint_readconfig

  !========================================================

  subroutine calc_bounds(lon, lat, lonb, latb)

    !> Calculates the boundaries between global grid-boxes. 
    !> Note that we assume that the boundaries lie half-way between the 
    !> points, both latitudinally and longitudinally, although 
    !> this isn't strictly true for a Gaussian grid.

    use glimmer_map_trans, only: loncorrect

    implicit none

    real(dp),dimension(:),intent(in) :: lon,lat    !> locations of global grid-points (degrees)
    real(dp),dimension(:),intent(out) :: lonb,latb !> boundaries of grid-boxes (degrees)

    real(dp) :: dlon

    integer :: nxg,nyg,i,j

    nxg = size(lon) ; nyg = size(lat)

    ! Latitudes first - we assume the boundaries of the first and 
    ! last boxes coincide with the poles. Not sure how to
    ! handle it if they don't...

    latb(1) = 90.d0
    latb(nyg+1) = -90.d0

    do j = 2,nyg
       latb(j) = lat(j-1) - (lat(j-1)-lat(j))/2.0
    enddo

    ! Longitudes

    if (lon(1) < lon(nxg)) then
       dlon = lon(1) - lon(nxg) + 360.d0
    else
       dlon = lon(1) - lon(nxg)
    endif
    lonb(1) = lon(nxg) + dlon/2.d0
    lonb(1) = loncorrect(lonb(1),0.d0)      

    lonb(nxg+1)=lonb(1)

    do i = 2,nxg
       if (lon(i) < lon(i-1)) then
          dlon = lon(i) - lon(i-1) + 360.d0
       else
          dlon = lon(i) - lon(i-1)
       endif
       lonb(i) = lon(i-1) + dlon/2.d0
       lonb(i) = loncorrect(lonb(i),0.d0)      
    enddo

  end subroutine calc_bounds

  !========================================================

  integer function check_mbts(timesteps)

    !> Checks to see that all mass-balance time-steps are
    !> the same. Flags a fatal error if not, else assigns that
    !> value to the output

    use glimmer_log

    implicit none

    integer,dimension(:) :: timesteps !> Array of mass-balance timsteps

    integer :: n,i

    n = size(timesteps)
    if (n==0) then
       check_mbts = 0
       return
    endif

    check_mbts = timesteps(1)

    do i = 2,n
       if (timesteps(i) /= check_mbts) then
          call write_log('All instances must have the same mass-balance and ice timesteps', &
               GM_FATAL,__FILE__,__LINE__)
       endif
    enddo

  end function check_mbts

  !========================================================

  subroutine check_init_args(orog_lats, orog_longs, orog_latb, orog_lonb)

    !> Checks which combination arguments have been supplied to
    !> define the global grid, and rejects unsuitable combinations

    use glimmer_log

    real(dp),dimension(:),optional,intent(in) :: orog_lats 
    real(dp),dimension(:),optional,intent(in) :: orog_longs 
    real(dp),dimension(:),optional,intent(in) :: orog_latb
    real(dp),dimension(:),optional,intent(in) :: orog_lonb 

    integer :: args
    integer,dimension(5) :: allowed=(/0,3,7,11,15/)

    args = 0

    if (present(orog_lats))  args = args + 1
    if (present(orog_longs)) args = args + 2
    if (present(orog_latb))  args = args + 4
    if (present(orog_lonb))  args = args + 8

    if (.not.any(args==allowed)) then
       call write_log('Unexpected combination of arguments to initialise_glint', &
            GM_FATAL,__FILE__,__LINE__)
    end if

  end subroutine check_init_args

  !========================================================

  subroutine check_input_fields(params, humid, lwdown, swdown, airpress, zonwind, merwind)

    use glimmer_log

    type(glint_params),              intent(inout) :: params   !> parameters for this run
    real(dp),dimension(:,:),optional,intent(in)    :: humid    !> Surface humidity (%)
    real(dp),dimension(:,:),optional,intent(in)    :: lwdown   !> Downwelling longwave (W/m$^2$)
    real(dp),dimension(:,:),optional,intent(in)    :: swdown   !> Downwelling shortwave (W/m$^2$)
    real(dp),dimension(:,:),optional,intent(in)    :: airpress !> surface air pressure (Pa)
    real(dp),dimension(:,:),optional,intent(in)    :: zonwind  !> Zonal component of the wind field (m/s)
    real(dp),dimension(:,:),optional,intent(in)    :: merwind  !> Meridional component of the wind field (m/s)

    if (params%enmabal) then
       if (.not.(present(humid)  .and. present(lwdown)  .and. &
                 present(swdown) .and. present(airpress).and. &
                 present(zonwind).and. present(merwind))) &
            call write_log('Necessary fields not supplied for Energy Balance Mass Balance model',GM_FATAL, &
            __FILE__,__LINE__)
    end if

    if (params%need_winds) then
       if (.not.(present(zonwind).and.present(merwind))) &
          call write_log('Need to supply zonal and meridional wind fields to GLINT',GM_FATAL, &
            __FILE__,__LINE__)
    end if

  end subroutine check_input_fields

  !========================================================

  subroutine accumulate_averages(params, temp, precip, zonwind, merwind, humid, lwdown, swdown, airpress)

    type(glint_params),              intent(inout) :: params   !> parameters for this run
    real(dp),dimension(:,:),         intent(in)    :: temp     !> Surface temperature field (celsius)
    real(dp),dimension(:,:),         intent(in)    :: precip   !> Precipitation rate        (mm/s)
    real(dp),dimension(:,:),optional,intent(in)    :: zonwind  !> Zonal component of the wind field (m/s)
    real(dp),dimension(:,:),optional,intent(in)    :: merwind  !> Meridional component of the wind field (m/s)
    real(dp),dimension(:,:),optional,intent(in)    :: humid    !> Surface humidity (%)
    real(dp),dimension(:,:),optional,intent(in)    :: lwdown   !> Downwelling longwave (W/m$^2$)
    real(dp),dimension(:,:),optional,intent(in)    :: swdown   !> Downwelling shortwave (W/m$^2$)
    real(dp),dimension(:,:),optional,intent(in)    :: airpress !> surface air pressure (Pa)

    params%g_av_temp    = params%g_av_temp    + temp
    params%g_av_precip  = params%g_av_precip  + precip

    if (params%need_winds) params%g_av_zonwind = params%g_av_zonwind + zonwind
    if (params%need_winds) params%g_av_merwind = params%g_av_merwind + merwind

    if (params%enmabal) then
       params%g_av_humid    = params%g_av_humid    + humid
       params%g_av_lwdown   = params%g_av_lwdown   + lwdown
       params%g_av_swdown   = params%g_av_swdown   + swdown
       params%g_av_airpress = params%g_av_airpress + airpress
    endif

    ! Ranges of temperature

    where (temp > params%g_max_temp) params%g_max_temp = temp
    where (temp < params%g_min_temp) params%g_min_temp = temp

  end subroutine accumulate_averages

  !========================================================

  subroutine accumulate_averages_gcm(params, qsmb, tsfc, topo)

    type(glint_params), intent(inout)   :: params     ! model parameters

    real(dp),dimension(:,:,0:),intent(in)  :: qsmb     ! flux of glacier ice (kg/m^2/s)
    real(dp),dimension(:,:,0:),intent(in)  :: tsfc     ! surface ground temperature (C)
    real(dp),dimension(:,:,0:),intent(in)  :: topo     ! surface elevation (m)

    integer :: nec    
    
    nec=params%g_grid%nec
    
    params%g_av_qsmb(:,:,0:nec) = params%g_av_qsmb(:,:,0:nec) + qsmb(:,:,0:nec)
    params%g_av_tsfc(:,:,0:nec) = params%g_av_tsfc(:,:,0:nec) + tsfc(:,:,0:nec)
    params%g_av_topo(:,:,0:nec) = params%g_av_topo(:,:,0:nec) + topo(:,:,0:nec)  
    !Topo accumulation and averaging is redundant, but retained for consistency with other input fields from GCM

  end subroutine accumulate_averages_gcm

  !========================================================

  subroutine calculate_averages(params)

    type(glint_params), intent(inout) :: params   !> parameters for this run

    params%g_av_temp    = params%g_av_temp   / real(params%av_steps,dp)
    params%g_av_precip  = params%g_av_precip / real(params%av_steps,dp)
    if (params%need_winds) params%g_av_zonwind = params%g_av_zonwind / real(params%av_steps,dp)
    if (params%need_winds) params%g_av_merwind = params%g_av_merwind / real(params%av_steps,dp)
    if (params%enmabal) then
       params%g_av_humid    = params%g_av_humid   /real(params%av_steps,dp)
       params%g_av_lwdown   = params%g_av_lwdown  /real(params%av_steps,dp)
       params%g_av_swdown   = params%g_av_swdown  /real(params%av_steps,dp)
       params%g_av_airpress = params%g_av_airpress/real(params%av_steps,dp)
    endif

  end subroutine calculate_averages

  !========================================================

  subroutine calculate_averages_gcm(params)

    type(glint_params),              intent(inout) :: params   !> parameters for this run
    
    integer :: nec
    
    nec=params%g_grid%nec    
    
    params%g_av_qsmb(:,:,0:nec) = params%g_av_qsmb(:,:,0:nec) / real(params%av_steps,dp)
    params%g_av_tsfc(:,:,0:nec) = params%g_av_tsfc(:,:,0:nec) / real(params%av_steps,dp)
    params%g_av_topo(:,:,0:nec) = params%g_av_topo(:,:,0:nec) / real(params%av_steps,dp)

  end subroutine calculate_averages_gcm

  !========================================================

  subroutine populate_output_flags(out_f,                     &
                                   orog_out,       albedo,    &
                                   ice_frac,       veg_frac,  &
                                   snowice_frac,   snowveg_frac,  &
                                   snow_depth,                &
                                   water_in,       water_out, &
                                   total_water_in, total_water_out, &
                                   ice_volume)

    type(output_flags),intent(inout) :: out_f

    real(dp),dimension(:,:),optional,intent(inout) :: orog_out        !> The fed-back, output orography (m)
    real(dp),dimension(:,:),optional,intent(inout) :: albedo          !> surface albedo
    real(dp),dimension(:,:),optional,intent(inout) :: ice_frac        !> grid-box ice-fraction
    real(dp),dimension(:,:),optional,intent(inout) :: veg_frac        !> grid-box veg-fraction
    real(dp),dimension(:,:),optional,intent(inout) :: snowice_frac    !> grid-box snow-covered ice fraction
    real(dp),dimension(:,:),optional,intent(inout) :: snowveg_frac    !> grid-box snow-covered veg fraction
    real(dp),dimension(:,:),optional,intent(inout) :: snow_depth      !> grid-box mean snow depth (m water equivalent)
    real(dp),dimension(:,:),optional,intent(inout) :: water_in        !> Input water flux          (mm)
    real(dp),dimension(:,:),optional,intent(inout) :: water_out       !> Output water flux         (mm)
    real(dp),               optional,intent(inout) :: total_water_in  !> Area-integrated water flux in (kg)
    real(dp),               optional,intent(inout) :: total_water_out !> Area-integrated water flux out (kg)
    real(dp),               optional,intent(inout) :: ice_volume      !> Total ice volume (m$^3$)


    out_f%orog         = present(orog_out)
    out_f%albedo       = present(albedo)
    out_f%ice_frac     = present(ice_frac)
    out_f%veg_frac     = present(veg_frac)
    out_f%snowice_frac = present(snowice_frac)
    out_f%snowveg_frac = present(snowveg_frac)
    out_f%snow_depth   = present(snow_depth)
    out_f%water_out    = present(water_out)
    out_f%water_in     = present(water_in)
    out_f%total_win    = present(total_water_in)
    out_f%total_wout   = present(total_water_out)
    out_f%ice_vol      = present(ice_volume)

  end subroutine populate_output_flags

  !========================================================

  subroutine compute_ice_sheet_grid_mask(ice_sheet_grid_mask, gfrac)

    !Calculate an array that contains the fractional area of ice+land.
    !This will ultimately be upscaled and sent to CLM, to indicate where
    !ice sheet can potentially exist.  For example, in the case of Greenland,
    !This mask is meant to define the contiguous island of Greenland, but not 
    !define ocean, or any other land points (e.g. Iceland) that fall within the
    !ice sheet grid, but are not represented by the ice sheet model.

    real(dp) ,dimension(:,:,0:),intent(in) :: gfrac               ! ice+bare land fractional area [0,1]
    real(dp) ,dimension(:,:), intent(out)  :: ice_sheet_grid_mask ! mask of ice sheet grid coverage


    !The following line sums gfrac over all exposed bare land (0-indexed) and ice
    !elevation bins.  Thus, it includes the contribution of bare exposed land to the
    !total fraction. 
    ice_sheet_grid_mask=sum(gfrac,3)

  end subroutine compute_ice_sheet_grid_mask

  !========================================================

  subroutine compute_icemask_coupled_fluxes(icemask_coupled_fluxes, ice_sheet_grid_mask, instance)

    ! Given an already-computed ice_sheet_grid_mask array, compute
    ! icemask_coupled_fluxes. The latter is similar to the former, but is 0 for icesheet
    ! instances that have zero_gcm_fluxes = .true. Thus, icemask_coupled_fluxes can be
    ! used to determine the regions of the world in which CISM is operating and
    ! potentially sending non-zero fluxes to the climate model.

    real(dp) ,dimension(:,:), intent(out)  :: icemask_coupled_fluxes ! mask of ice sheet grid coverage where we are potentially sending non-zero fluxes
    real(dp), dimension(:,:), intent(in)   :: ice_sheet_grid_mask    ! mask of ice sheet grid coverage
    type(glint_instance), intent(in)       :: instance               ! the model instance


    if (instance%zero_gcm_fluxes == ZERO_GCM_FLUXES_TRUE) then
       icemask_coupled_fluxes(:,:) = 0.d0
    else
       icemask_coupled_fluxes(:,:) = ice_sheet_grid_mask(:,:)
    end if

  end subroutine compute_icemask_coupled_fluxes

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

end module glint_main

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
